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We investigate the net-baryon number fluctuations across the chiral phase transition 
at finite density in the strong coupling and chiral limit. Mesonic field fluctuations are 
taken into account by using the auxiliary field Monte-Carlo method. We find that the 
higher-order cumulant ratios, Sa and kct^, show oscillatory behavior around the phase 
boundary at /r/T > 0.2, and there exists the region where the higher-order cumulant 
ratios are negative. The negative region of is found to shrink with increasing lattice 
size. This behavior agrees with the expectations from the scaling analysis. 
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1. Introduction 


Fluctuations of conserved charges are promising observables in relativistic heavy-ion colli¬ 
sions in search for QCD phase transition . In the first phase of the Beam Energy Scan 

(BES I) program at Relativistic Heavy Ion Collider (RHIC), net-proton Q], as a proxy of 
net-baryon 0,0, and net-electric charge event-by-event fluctuations 0 have been measured 
in Au-|-Au collisions for a broad energy range from ^snn = 7.7 to 200 GeV. The success of 
the statistical thermal model 0 indicates the event-by-event fluctuations of the multiplicity 
of conserved charges can be regarded as particle number fluctuations in the grand canonical 
ensemble specified by temperature, volume and chemical potential. 0 
The particle number fluctuations of the conserved charge, i.e., susceptibilities or cumulants, 
reflect the property of the phase in QCD lOl], which is expected to exhibit a rich structure 


ll|. Especially, if QCD has a critical point (CP) 1^, the susceptibilities and higher-order 


cumulants diverge at the critical point owing to the divergent correlation length 13|,ll^, thus 
one expects anomalously large fluctuations would be observed if the system passes around 
CP. In heavy-ion collisions, the correlation length cannot diverge due to finite size of the 
system, but these fluctuations are expected to remain sensitive to the remnant criticality 
of the system around CP 0. Since this sensitivity implies the tail of the event-by-event 
multiplicity distribution has importance in the behavior of the cumulants la], the analysis 
becomes statistically demanded. 

The measured higher-order cumulant ratios of net-proton number, Sa = Xp IXp 

= Xp^/x}p \ show non-monotonic behavior as a function of the incident energy, where 
(T^, S and K are referred to as the variance, skewness and kurtosis, respectively. At around 
= 20 GeV, the data show < 1 below the expected value in the Skellam distribution 


The decrease of in the net-proton number can be the signal of the critical behavior. 
According to the theoretical arguments 16|,ll7|, of the net-baryon number can be reduced 
by critical behavior from universality. In QCD, the expected universality class depends on the 
quark masses. Eor two-flavor massless quarks with axial anomaly, the QCD phase transition 
at finite T (// = 0) belongs to the universality class of 3d 0(4) symmetric spin model [H, 
0 , and the second order transition at low /i may turn into the hrst order transition at 
the tricritical point (TCP). At physical quark masses, finite T phase transition would be 


crossover 


20|, and the fluctuations are governed by the approximate chiral symmetry. The 


pseudo-critical line at low /i may be connected with the hrst order transition line at CP, 
whose criticality is expected to belong to Z(2) universality class 


12 


mM 


While the universality argument dictates the singular behavior of the thermodynamic 
quantities close to the phase transition, the actual magnitude of the huctuations are smeared 
by the hnite volume effect in addition to the hnite quark mass. Finite volume makes the 
transition smoother. We need to perform hnite size scaling analysis to observe precise critical 
behavior 20|. The critical behavior depends on the relative strength of the singular part to 


the regular (non-singular) part of the free energy density (or its derivatives). The regular part 


^ This argument holds for a limited acceptance. If one covers entire phase space in heavy ion 
collisions, those fluctuations is of course absent. One^an invoke more informations by changing the 
rapidity acceptance beyond the equilibrium regime [9|. 
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would also mask the expected singular behavior from criticality 16|]. Thus, explicit calcula¬ 
tions are desirable to see how much criticality can be present in the fluctuation observables, 
and to pin down the origin of the observed decrease of Z(2) critical behavior around 
CP 


14l |. the remnant of the 0(4) chiral phase transition 1^, or other mechanisms. 


Lattice QCD (LQCD) is the most powerful non-perturbative approach to QCD based 
on Monte-Carlo (MC) simulations, and is successful at vanishing or low baryon chemical 
potential. For instance, higher-order cumulants have been studied at zero chemical potential 


22 , 


23 , [24] and at finite chemical potential 23, [ 23 , |2a, [ 23 , l28|. At large baryon chemical 


potential, however, LQCD faces the notorious sign problem which makes it difficult to carry 
out MC simulations owing to complex fermion determinant. There are many attempts to 
circumvent the sign problem 


29, 


3C, 


31 


32 I . I 33 I . I 34 I . I 35 I . l36l | . Most of the methods evading 


the sign problem are reliable only in limited circumstances such as p/T < 1, small volume, 
or heavier quark masses than physical one, depending on the method. Conclusive results at 
physical point have not been obtained yet. 

Owing to the sign problem in LQCD, most of fluctuation studies at nonzero density have 
been carried out by using chiral models [la, [la, l37[, [33, l39[, |40[, |4]|. It has been pointed out 
that taking into account fluctuations around the phase transition is essential to correctly 
describe the behavior of the fluctuation of conserved charges 15l.ll6l.l37l|. In chiral models, this 
has been done by implementing functional renormalization group (FRG) 43, [4^. This is a 
natural consequence of the fact that the critical behavior of the conserved charge fluctuations 
is governed by the critical exponent of the specific heat a which requires beyond mean-field 


treatment 


la, y,|37j,|38[,|39(|. 


One of the possible alternative ways to attack the finite density re gion is the strong cou¬ 
pling approach of LQCD with fermions 43, lia, |46|, |47|, |48|, |49|, [soj, Isil, [53, 


5S, 


59, 6C, 


53[, [53, |55[, |56[, [53, 

6ll | . The strong coupling approach with fermions could reduce the severity of the 


sign problem and is applied to investigate the chiral phase transition. Recently, theoretical 
frameworks including flnctnation effects in the strong coupling limit (SCL) have been devel¬ 
oped: auxiliary field Monte-Carlo (AFMC) method 4^ and the monomer-dimer-polymer 
(MDP) simulation 45, 46, 47, 3,1^. These two methods give consistent phase boundaries 
of the chiral phase transition in the chiral limit. Althongh SCL is the coarse lattice limit, we 
might expect that the observables related to the phase transition is not so sensitive to the 
coarse lattice spacing since long-wave dynamics dominates the property of the phase tran¬ 
sition, inclnding the influences on the higher-order cumulants |39l ]. Fnrthermore, one can 
take the chiral limit in which the chiral transition becomes second order. Thus, the singular 
behavior associated with the phase transition is smeared by the finite size effect only. As a 
result, the divergent part of the higher-order cumulants, which can be described by relevant 
scaling fnnction of the singnlar part of the free energy density , conld be replaced by sign 
changes across the transition (la . HQ- 
In this article, we focus on the fluctuations of net-baryon number and discuss the higher- 
order cumulant ratios in the strong coupling limit of LQCD by utilizing the AFMC method 
to take the fluctuation effects into account. We here consider the LQCD action with one 
species of unrooted staggered fermion in the chiral limit. We demonstrate that the higher- 
order cnmulant ratios, Sa and kct^, show oscillatory behavior, and there exists the region 
where higher-order cnmnlants are negative in the strong conpling and chiral limit on not 
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too small lattices at fi/T > 0.2. We also discuss the lattice size dependence of the negative 
region. 

This paper is organized as follows. We briefly provide the formalism in Sec. [2j Onr main 
results on the higher-order fluctuations of the net-baryon nnmber will be presented in Sec. [31 
We summarize our paper in Sec. 01 Some formnlae for the net-baryon number cumulants 
used in Sec. El are fonnd in the appendix. 


2. Lattice QCD in the strong conpling limit with fluctuations 

We consider a lattice QCD action with one species of unrooted staggered fermion in the 
d{= 3) -|- 1 dimensional anisotropic Enclidean spacetime with W temporal and L spatial 
lattice sizes. This LQCD action has remnant chiral symmetry U{1)r x U{1)l in the chiral 
limit (mo —)• 0). We set the nnmber of the color Nc = 3 and the lattice nnit a = 1 thr ong hont 
this paper. In the following, we briefly summarize the formalism developed in Ref. 441. 


2.1. Effective action in the strong coupling limit with fluctuations 

In the strong conpling limit (SCL) g oo, we conld ignore the plaquette terms, which are 

proportional to 1/g^, so the partition function and the action of the lattice QCD become 

^SCL = [ V[x, X, Uu] , 


SF=i:^[Vff-V- 


1 


d 


X j = l 


( 1 ) 


XxUj,^X^^j - X^+jUl^Xx , (2) 


=7e^''-^^'^^Xxlto,xX:,+o , Vx =je ^''■^^'^^Xa^+d^lxXx , = XxXx , (3) 

where Xx and represent the staggered quark field and the link variable, respectively, 
and Vff and Mx are mesonic composites. The staggered sign factor pj^x = (— 
is related to the gamma matrices in the continunm limit. The quark mass mg is taken to 
be zero throughout this paper. Quark chemical potential g is introduced together with the 
physical lattice spacing ratio /(y) = where 7 is the anisotropy parameter. We 


adopt /(y) = y^ following the arguments in Refs. 44l. l48l . l54l . l55 1 


We obtain the effective action of the auxiliary fields, , after the followiM three steps 


First, by integrating ont spatial link variables [^, [Sl], j^, [SSj, [bj, [ba, l56l. l59l. [GOj in the 
leading order of the strong conpling and 1 /d (large-dimensional) expansion ^ , one finds a 
convenient expression for the effective action. 


Sen = ^ E ['^ 


- U- - 


1 

4W 




x+j 


+ mo'^Mx 


(4) 


X,J 


It should be noted that spatial baryonic hopping terms are not included in Eq. (0|) as they 
are higher order in the 1/d expansion. By comparison, we exactly integrate out temporal 
link variables later, then the temporal baryonic hopping effects are taken into acconnt. 

In the second step, we transform the effective action to the fermion-bilinear form by using 
the extended Hnbbard-Stratonovich (eHS) transformation 59], [60|, 

e»AB ^ f ^^^^^-a{l^-[A+B)/2]^ + [<l>-i(A-B)/2r}+aAB 




( 5 ) 
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where 'll) = (p + icj) and dij) dtp* = dUeip dYmip = dipdp. The four-Fermi interaction terms in 
Eq. (jH) are diagonal in the momentum space, and are separated into two parts based on the 
momentum regions: the positive modes (/(k) = > C)) and the negative (/(k) < 

0 ) modes. 

L3 


■ E = -^ E /(X) " 1 -. 

^ k,T 


L3 

'4Nr 


E 


/(k)(Mk,,M_k.. - Mk,,M_k,, 


( 6 ) 


k,r,/(k)>0 


where = Xlk ^-^k.D k = k-|- (vr, 7 r, 7 r), and /(k) = —/(k). The effective action 

after the eHS transformation of Eq. ([ 6 ]) reads 


qEHS _ ^ 
-^eff “2 

X 


[Vp~ - Ej. ] + ^ nixMx: + 


4Nr_ 


E /(k) 


Tk,'; 


|2 I I |2 

I T Pk^rl 


k,'r,/(k)>0 


rrix =mo + 


— E 

j 


where cJa, = Z]k, 


/(k)>0 ' 


a + zevr) 

:*^'^CJk,r and TTx = X]k, 


X+] + 


(7) 

( 8 ) 


,/(k)>o(“l)’'®*’^ '^^k,r- i^Ote riiar 0-_k,r - ‘^k,r 
are bosonized at a time. At small k, they are 
also regarded as the chiral (cr) fields and the Nambu-Goldstone (vr) fields, respectively, since 
the staggered fermion identifies the spin and flavor by specifying space-time position. The 
sign factor, Sx = l—\l^o+xi+x 2 +x 3 ^ corresponds to 75 ( 81*75 in the spinor-taste space in the 
continuum limit. 

In the third step, by integrating over the Grassmann and temporal link {Uq) variables 


Note that (T_k ^ = cr,* 


and vr_k,T = i^kr> ®i^ce ±k terms in Eq. 


IM, l55( , the partition function and the effective action are reduced to 

Zay = J T>[crk,r,7rk,r] , 

EV(k) 


cAF _ 
<^65 — 


E 

k,r,/(k)>0 


4iV, 


Fk,-i 


|2 I 1.^ |2 


^logi?(x) , 


R{x.) =AAr^(x)^ - 2AAr^(x) + 2cosh{Ncn/T) 


(9) 

( 10 ) 

( 11 ) 


where P [cTk,r, 7rk,r] = nk,r,/(k)> 0 '^^k,rdo-k^7%; 


dvrj)^ and T = In the last line, we 

use a recursion formula to obtain [^,1^,1^. In the cases where 'mx=[^^r) is static, 

we And = 2 cosh(AT- arcsinh (nix/'y))- It should be noted that the action represents the 
confinement phase, since the baryonic chemical potential appears as cosh(Ac/r/T), which can 
be understood as baryonic contribution 6 II], similarly to the PNJL model with vanishing 
Polyakov loop 621. 

We can now perform the Monte-Garlo integral over the auxiliary fields ((7k,T) ^Ck.r) using the 
effective action 5^^. We generate Monte-Garlo configurations based on the phase quenched 
action at finite p and T and calculate observables by the reweighting method; 

{O exp(i6»))pq 


{ 0 ) = 


( 12 ) 


(exp(i6'))pq 

where 9 = —Im(5^^) and (• • • )pq denotes the phase quenched average. Through this AEMC 
method, we can numerically take account of fluctuation effects. 
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It should be noted that we have a sign problem that stems from the bosonization procedure, 


but it is not severe and we can investigate the QCD phase diagram as done in Ref. 4^ 


The reason for the milder sign problem may be understood as follows. First, the auxiliary 
field effective action is obtained by integrating out the link variables, then it contains only 
the color singlet field. Color singlet states are, in general, closer to the energy eigenstates 
than colored states, and we expect smaller phases in the path integral. Next, there is no 
sign problem in the mean field approximation, where Cx is assumed to be constant and tTx 
fields are neglected, then the sign problem is not severe with small |fc|. In the case where 
the long wave physics dominates, only small |k| auxiliary fields are relevant, thus the vr 
field can be regarded as almost constant. As seen in Eq. ([8]), the complex phase of one site 
has an opposite sign to that of the nearest neighbor sites. Therefore, we could expect that 
the complex phase on one site coming from the bosonization is canceled out by the nearest 
neighbor site contributions as long as we study the long wave phenomena. As a result, we 
have a milder sign problem and can directly generate MC configurations at finite and T 
as long as the lattice size is not very large. In actual calculations, high |k| modes are not 
negligible, and the average phase factor (exp(i0))pq is suppressed; @ (exp(i0))pq > 0.9 and 0.4 
for ^l/T < 0.8 on the 4^ X 4 and 8^ x 8 lattices around the phase boundary, respectively 44l |. 

In Fig. [H we show the chiral condensate (u = ^i«=o,t/^t)) and the logarithm of 

the baryon number susceptibility = d'^\ogZ/d{NcfJ,/T)‘^/{VT^)) in the chiral limit 
on a 6^ X 6 lattice as a function of T/Tc, where Tc(~ 1.46862) denotes critical temper¬ 
ature at ^ = 0 on a 6^ X 6 lattice defined by the peak position of chiral susceptibility 
iXcr = log Z/drriQ/{L^Nx)) j^. We can see the chiral phase transition behavior from the 
decrease of the chiral condensate with increasing T. The decrease of the chiral condensate 
becomes steeper with increasing /i/T. The baryon number susceptibility has a sharp peak at 
high /i/r. These behaviors suggest the influence of the singularity at the tricritical point. It 
should be noted that the baryon number susceptibility at zero chemical potential decreases 
at high T. This trend is not consistent with standard lattice QCD results and would be an 
artifact of the truncation of the action in Eq. (SD, so we only focus on the vicinity of the 
phase transition. 


2.2. Some remarks on the universality class 

The staggered fermion has a remnant chiral symmetry in the chiral limit and has 0(2) sym¬ 
metry as long as the lattice spacing is coarse. We introduce the auxiliary fields to maintain 
the original symmetry, 0(2). The auxiliary field action in Eq. (jlOp is invariant under 
the chiral transformation. 


(^k\ (o-'i\ ^ /cosT -sinT\ / ak\ 

yk) l^sinT cosT J J 

which corresponds to the transformation of the staggered fermion fields 


Xx 


x'x = 


Xx 


x'x = 


(13) 


(14) 


^ The average phase factor indicates the severity of weight cancellation. We have no weight 
cancellation when (exp(i0))pq = 1. 
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chiral condensate, 6 X6 



T/T, 


Fig. 1 Chiral condensate and logarithmic 
6^ X 6 lattice for various ^/T lines. The cross 
diamond, and pentagon denote /r/T = 0.0,0.1 


3 

Baryon number susceptibility, 6 X6 



T/T, 


scale of baryon number susceptibility on a 
mark, square, circle, triangle, lower triangle, 
0.2,0.3, 0.5,0.6 and 0.8, respectively. 


In the framework of the MDP in the strong coupling limit, the critical exponent was studied 
and that values are consistent with 0(2) [46[. It should be noted that we need careful 
discussions about the difference in symmetry since the critical exponents are different in 
0(2), 0(4), and Z(2). We will mention the relation between the 0{N) scaling function and 
critical behavior in Sec. 


3. Net-baryon number fluctuations in the strong coupiing iimit 

3.1. Net-baryon number cumulants 

In this section, we present results on the higher-order cumulant ratios of the net-baryon 
number. The n-th order cumulant of the net-baryon number in the grand canonical ensemble 
is given by 


yW = 


1 5” log Z 


A = Ncy/T . 


(15) 


FT3 dfi^ 

To discuss the effects of the phase transition on the net-baryon number fluctuations, it is 
convenient to take the ratio of a higher-order cumulant to the second order one 63, since 
the volume factor, V, in Eq. (jl5|l cancels. We show the results of the following cumulant 
ratios, the normalized skewness and kurtosis. 


= = (16) 


The skewness S and kurtosis k probe the asymmetry with respect to the mean and peaked- 

( 2 ) 

ness of the underlying probability distribution, respectively. The normalization by Xij. 
implies that one sets a reference, since = 1 for the Skellam distribution which describes 
the distribution of the difference ni — n 2 , where ni and n 2 follow the Poisson distribution. 
Then it corresponds to the free hadron gas with Boltzmann approximation. 

The detailed method for the calculation of the cumulants are summarized in Sec. As 
discussed in Sec. El observables have an imaginary part due to the high |k| modes of vr^^r, 
so we take the real part of the observables and error bars in this article. The absolute mean 
value of maximum imaginary part for the normalized kurtosis and skewness are about 70 
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and 2 for /i/T = 0.8 and 0.5 and 0.04 for /r/T = 0.2, respectively, which are smaller than 
the peak heights and valley depth of the real part as seen in Figs. [3] and 01 It should be 
noted that the high temperature behavior of the cumulants would bare the artifacts of the 
present treatment. Figures 0] and 0] indicate that Sa stays negative in the large chemical 
potential region and can also take small negative values depending on ^/T at high T. 
For a free ideal bary on gas we expect Sa —)■ + vr^) and ^ + vr^) at high 


temperature 2^, l63l . [65| . These differences at high temperature may be due to the truncation 


in the 1/d expansion. Spatial baryon hopping term is missing in the leading order in the 
1/d expansion, then the fermion momentum integral generates only a volume factor. As a 
result, pressure is proportional to T rather than in the Stefan-Boltzmann behavior at 
high temperature. Since we are interested in the transition region, where we expect that the 
long wave mesonic fluctuations dominate, we focus on the cumulant ratios around the phase 
boundary in the later discussion. 


3.2. T and ^ dependence 


Sa, 6^ X 6 Ko^, 6^ X 6 




Fig. 2 Normalized skewness (left panel) and kurtosis (right panel) on the various p/T lines 
on a 6^ X 6 lattice The cross mark, circle, triangle, lower triangle, diamond, and pentagon 
indicate p/T = 0.0, 0.2, 0.3, 0.5, 0.6, and 0.8 line, respectively. As for the skewness on 
the p/T = 0.8 line, we rescale results by factor 3. For the kurtosis on the p/T = 0.5,0.6,0.8 
lines, we rescale results by factor 2, 3, 20, respectively. 


In Fig. [21 we show the results of Sa and on a 6^ lattice at several p/T as a function of 
T/Tc. One sees the strong dependence of these higher-order cumulants on both temperature 
and chemical potential. In the low chemical potential region, both Sa and are positive 
and have a broad peak. 

As p/T increases, characteristic structures emerge. Sa increases with T, then exhibits a 
strong positive peak followed by sudden decrease to large negative value. Further increase of 
temperature leads to a constant small value. The negative region of the skewness, indicating 
the critical behavior 4^, starts to appear between p/T = 0.3 and p/T = 0.5. shows 


similarly moderate increase with temperature and a sharp positive peak. Then it turns to 
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Sa, [irT=Q.2 


So, |,i/T=0.8 




Fig. 3 Normalized skewness on the /r/T = 0.2 (left panel) and /r/T = 0.8 (right panel) 
line on 4^ X 4 (the cross mark), 6^ x 4 (the diamond mark), and 6^ x 6 (the circle mark) 
lattices. A horizontal line around zero denotes the mean filed value at high temperature 
(see text). {Sa ~ —0.0220 for /r/T = 0.2 and Sa ~ —0.5442 for /r/T = 0.8.) We can find the 
negative region in the high chemical potential region. 


be negative with a sharp valley but becomes positive again with a somewhat milder peak. 
Such an oscillatory behavior appears in ///T > 0.2. Both cumulant ratios have one negative 
valley, and the skewness (kurtosis) has one (two) positive peak(s) at high chemical potential. 
Comparing with the phase boundary 4^, one Ends these behaviors appear around the phase 
boundary (see below). In the following subsections, we discuss the behaviors in more details 
by looking at lattice size dependence. 


3.3. Lattice size dependence 

In Figs. [3] and 01 we show the results of the normalized skewness Sa and the normalized 
kurtosis respectively, on the fixed p/T lines, fi/T = 0.2 (left panel) and p,/T = 0.8 
(right panel), on 4^,6^ x 4 and 6^ lattices. Horizontal lines show the mean field values at 
high temperature, where the chiral condensate vanishes. We will use these mean field values 
in Sec. 13.51 to remove the artifact as discussed in Sec. 13.11 

Lattice size dependence of Sa is moderate at low p/T and prominent at large p/T. From 
Fig. 01 one sees that the qualitative structure of the temperature dependence of Sa does 
not strongly depend on the lattice size. For p/T = 0.2, the skewness is positive around the 
phase transition region and has a peak. While the peak height does not show statistically 
signihcant dependence on the lattice size, one finds the temperature dependence becomes 
slightly stronger in the largest one, 6“^. The lattice size dependence becomes prominent at 
large chemical potential p/T = 0.8. The skewness has one positive peak and one negative 
valley at p/T = 0.8. The widths of the peak and the valley become narrower on a 6^ lattice. 

Similar lattice size dependence is found for na"^. In the low chemical potential region 
[p/T = 0.2), the smallest lattice size 4^ does not exhibit the characteristic structure seen in 
larger lattice cases, but one sees a moderate decrease around T/Tc ~ 0.9 following a peak 
at lower T. The broad negative valley appears on 6^ x 4 and 6^ lattices, implying that one 
should not take too small lattice size to see the influence of the critical fluctuations. We hnd 
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Ka^, |a,/T=0.2 ko^, ^iyT=0.8 



T/T^ T/T^ 


Fig. 4 Normalized kurtosis on the ^/T = 0.2 (left panel) and ^/T = 0.8 (right panel) line 
on 4^ X 4 (the cross mark), 6^ x 4 (the diamond mark), and 6^ x 6 (the circle mark) lattices. 
The thin horizontal line around zero denotes the mean field value at high temperature (see 
text), (kct^ ~ —0.1033 for /r/T = 0.2 and ~ —0.0251 for ^/T = 0.8.) 


that the oscillatory behavior with two positive peaks and one negative valley and its lattice 
size dependence becomes prominent for ///T = 0.8. Again, as lattice size becomes larger, 
the peak height tends to be larger value and the width of the oscillatory region becomes 
narrower. 


3.4- Relation to 0{N) scaling 

The behavior of the higher-order net-baryon number cumulants around chiral phase 
transition are expected to be described by the property of the 0{N) scaling function 


Ifil . IBOI . 1671 . 1681 . [69| . According to Ref. [161] , the first divergent cumulant at // = 0 in the 


thermodynamic limit is the sixth order one for N = 2 and 4. This result comes from the 
fact that the critical exponent of th e sp ecific heat a is negative in these cases, a ~ —0.0147 
in 0(2) 6^ and a ~ —0.21 in 0(4) 67|. For positive a, the divergence appears at the fourth 
order cumulants while it has a cusp for a < 0 15|, l37(| . 

For a nonzero /x, the first divergent cumulant is the third order one thus Sa diverges in 
the thermodynamic limit for 3d 0(2) and 0(4) universality class. Once the singular part is 
smeared by the explicit breaking (nonzero current quark mass) or finite volume effects El, the 
leading singular contribution is suppressed by the small multiplicative factor —(2xt^)(u/T)’^ 
where Kq denotes the curvature of the phase boundary near // = 0 in the chiral limit |l6l |. 
Then, whether one sees the effects of critical fluctuations in the higher-order cumulants or 
not depends on the value of /x/T and magnitude of the regular part of the free energy density 
and its derivatives. For instance, Sa positively (negatively) diverges in the thermodynamic 


^ Recalling the behavior of the order parameters becomes moderate not only by finite mass but 
also finite size effects, we could guess that the singular part of the cumulants is also masked by the 
finite size effects. By using models, we can find the characteristic behavior of physical quantities 
becomes smoother when the system is finite (70|. In the framework of 3d Ising model, finite size 
results are studied in Ref. (tH and the peak heights of Sa and na^ increase with larger volume due 
to the correlation length T The relations between higher-order net-bax^yon number cumulants and 
the correlation length in 3d 0(4) spin model are pointed out in Ref. [T^. 
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limit when approached from below (above) the phase transition. The negative divergent 
contribution is replaced by a large negative value in the presence of the explicit breaking 
or hnite volume effects. Then, the appearance of the negative Sa depends on whether the 
singular part overcomes the regular part. Similar argument applies to and the absence 
of the negative region in in the smallest lattice (Fig. [1]) is a consequence of the finite 
volume effects. Since exhibits the oscillation around the phase transition along constant 
/r/r including the narrow negative region between the two positive peaks, one expects 
positively diverges both from below and above the phase transition and negative region 
would disappear in the thermodynamic limit. The lattice size dependence of the negative 
region meets this expectation. This is in contrast to the case of the finite quark mass. 
With finite quark mass, one finds the negative kurtosis region even in the thermodynamic 
limit as a remnant of the divergence in the chiral limit la, l37 1. 


3.5. Kurtosis in the QCD phase diagram in the strong coupling limit 



Fig. 5 The negative kurtosis region in the QCD phase diagram in the chiral and strong 
coupling limit. The shaded area denotes the region where the is both negative and 
smaller than the mean field value on a 6^ lattice (more detail explanations are in the text 
and footnote). The vertical axis is the temperature parametrized as T = and the 

horizontal axis is the quark chemical potential. Both T and g are normalized by using Tc, 
the critical temperature at /z = 0 on a 6^ lattice. The phase boundaries on 4^ X 4,6^ X 4, 
and 6^ X 6 lattices are determined by the peak position of the chiral susceptibility {xa = 
(9^ log Z/(9mQ/(L^W)) for g/T < 0.8 and the effective potential analysis dehned by Fgff = 
(5'efr) /{KtL^) for g/T > 1.0 4^. A gap is the expected boundary between the would-be 1st 
order and the would-be 2nd order transition. 


In Fig. [5l we show the negative kurtosis region in the chiral limit (mo —)• 0) on a 6^ lattice. 
We also show the phase boundary determined by the peak position of the chiral susceptibility 
{X(T = d'^ log Z/dmQ/{L^N t)) for g/T < 0.8 and by the cross point of the effective potential 
dehned by Fefj = (5eff) /{Nt-L^) for g/T > 1.0. The transition would be the second order for 
g/T < 0.8. While the numerical hnite size scaling analysis cannot rule out crossover due to 
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the statistics 44], we expect the second order transition from symmetry arguments [I3|. We 
guess that the tricritical point exists at /i/T > 0.8. 

We here define the negative kurtosis region where is smaller than the mean field results 
at high T, in order to reduce the artifact at high T as discussed in Sec. 13.11 It should be 
noted that we do not take account of error bars to define the region here. The negative 
kurtosis region appears from about n/T = 0.2, and the region is the largest at ^/T = 0.30. 
By comparison, the kurtosis is positive at /r/T < 0.2. This might be due to the suppression 
of the singular contribution by the factor (/i/T)” and the regular part dominates over the 
singular part la]. 

The negative region almost coincides with the phase boundary. As discussed in Secs. 13.31 
and l3.4l the region will shrink to the phase boundary and finally vanish in the thermodynamic 
limit. Although the dependence of the scaling function is different between the finite size 
effect and the symmetry breaking term, one could expect that the negative region also 
appears with finite mass in the strong coupling limit. 

We conclude that we find the negative skewness or kurtosis region in the AFMC method 
in the chiral and strong coupling limit due to the finite volume effects as for the kurtosis, as 
a consequence of the critical fluctuations around the phase boundary incorporated through 
the AFMC method. 


4. Summary 

We have investigated the higher-order cumulant ratios of the net-baryon number, Sa = 

ill the strong coupling {g —)• oo) and the chiral limit (mo —)• 0) 
of QCD in the leading order of large dimensional expansions. Mesonic fluctuation effects 
are taken into account by making use of the auxiliary field Monte-Carlo (AFMC) method. 
We find the negative skewness and the kurtosis region around the boundary of the chiral 
phase transition. The skewness and kurtosis exhibit characteristic temperature dependences 
influenced by the critical fluctuations. 

The skewness is found to be negative near the phase boundary for large ^i/T. The oscilla¬ 
tory behavior around the phase boundary seems to be consistent with the expectations from 
the finite volume effect such that the positive (negative) peak at lower (higher) temperature 
diverges in the thermodynamic limit [l^. Similarly, the kurtosis has one negative valley 
between two positive peaks for large /r/T. With increasing lattice size, the negative valley 
is found to shrink, as anticipated for the finite volume effect. Thus, we expect two positive 
peaks diverge and the negative region disappears in the thermodynamic limit 16 j. 

One of the important next steps could be to study the effect of finite mass. This can be 
carried ou t by applying the AFMC method. Another important step is to see the finite size 


effects 20l. 1 701.1731] and determine the negative kurtosis area with or without bare quark mass. 
Whether one could have such region in lattice QCD is a mandatory question. As discussed 
above, it will depend on whether the smeared singular contribution by finite volume and/or 
explicit symmetry breaking term can overwhelm the regular contribution 16|, l70(], which 
corresponds to hadron resonance gas in finite temperature QCD below the phase transition. 


Along the g/T = 0.3 line at high T, the mean field value and AFMC results are so close to 
each other and it seems that the AFMC results oscillate around the MF result. Thus we define the 
boundary of the negative kurtosis region as the first intersection of the AFMC and MF results. 
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Since our present formulation ignores the spatial baryon hopping, we would expect that the 
regular contribution might be smaller than the realistic case thus our results could serve as 
a lower limit for the value of chemical potential where < 0. 
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A. Higher order derivatives of the baryon chemical potential 

In this appendix, we show the higher order derivatives with respect to dimensionless chemical 
potential fi{= Ncn/T). In the following, the action S and the partition function Z are denoted 
as and Z^-p, respectively. Then, higher order derivatives are given as 
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where fi = Nc^i/T. The derivatives of the action with respect to fi are given as 
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f{k) [|(Tfc,rP + kfc,rp] - Ea, ^og 7?(x) and i?(x) = X%{x) - 2 Xn{x) + 


where S = Efc,r,/>o Tv: 

2 cosh p. 

In this analysis, we have a complex phase coming from the fermion determinant, so the 
observables also take a complex value even if the observables are real in essence. Invoking the 
cancellation of the imaginary part in the case of higher statistics in principle, we can take 
the real part of the observables. We here use the jackknife method in order to evaluate the 
statistical errors 7^. Since the observables have an imaginary part, we evaluate the error 


with complex phase and take the correlation between complex observables and the complex 
phase into account. Then, we take the real part of the obtained error bars. 

When we evaluate the fourth derivative of logZ with respect to p, Xf^\ (the function 
form is like (O) = (a) + (6) (c) + (d)^), we generate jackknife samples of each expecta¬ 
tion value ((a)bin 5 (^)birn ('^)birn (h)bin)' Then, we calculate the jackknife samples of the 
observable ((Cl)bin = (®)bin + (^)bin ('^)bin + (h)bin)' Finally we calculate expectation value 
and error bars of the observable by using (Cl)bin- When we calculate the cumulant ratios 
(k(T^ = Sa = we use the same technique as well. 
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